﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// The Normal distribution is a distribution over the real numbers. The distribution
    /// is parameterized by a mean and standard deviation. The pdf is
    /// 
    ///     p(x) = \frac{1}{\sqrt{2 \pi} stddev} \exp{-0.5 (x - \mu)^2 / stddev^2}.
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class Normal : IContinuousDistribution
    {
        // The mean of the normal distribution.
        private readonly double mMean;

        // The standard deviation of the normal distribution.
        private readonly double mStdDev;

        /// <summary>
        /// Constructs a standard normal (mean = 0, standard deviation = 1).
        /// </summary>
        public Normal() : this(0.0, 1.0)
        {
        }

        /// <summary>
        /// Constructs a normal distribution.
        /// </summary>
        /// <param name="mean">The mean for the normal distribution.</param>
        /// <param name="stddev">The standard deviation for the normal distribution.</param>
        /// <exception cref="ArgumentException">If <paramref name="stddev"/> is negative.</exception>
        public Normal(double mean, double stddev)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(stddev);
            }
            
            mMean = mean;
            mStdDev = stddev;
            RandomNumberGenerator = new System.Random();
        }

        #region IDistribution Members

        /// <summary>
        /// Gets the mean of the distribution.
        /// </summary>
        /// <value>The mean of the distribution.</value>
        public double Mean
        {
            get { return mMean; }
        }

        /// <summary>
        /// Gets the standard deviation of the distribution.
        /// </summary>
        /// <value>The standard deviation of the distribution.</value>
        public double StdDev
        {
            get { return mStdDev; }
        }

        /// <summary>
        /// Gets the variance of the distribution.
        /// </summary>
        /// <value>The variance of the distribution.</value>
        public double Variance
        {
            get { return mStdDev*mStdDev; }
        }

        /// <summary>
        /// The entropy of the distribution.
        /// </summary>
        /// <value></value>
        public double Entropy
        {
            get { return System.Math.Log(mStdDev * Math.Constants.Sqrt2PiE); }
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }

        #endregion

        #region IContinuousDistribution Members

        /// <summary>
        /// The mode of the distribution.
        /// </summary>
        /// <value></value>
        public double Mode
        {
            get { return mMean; }
        }

        /// <summary>
        /// The median of the distribution.
        /// </summary>
        /// <value></value>
        public double Median
        {
            get { return mMean; }
        }

        /// <summary>
        /// Evaluates the probability density function for a the normal distribution.
        /// </summary>
        public double Density(double x)
        {
            double delta = (x - mMean) / mStdDev;
            return System.Math.Exp(-0.5 * delta * delta) / (Math.Constants.Sqrt2Pi * mStdDev);
        }

        /// <summary>
        /// Samples draws a random sample for the distribution.
        /// </summary>
        /// <returns>A random number from this distribution.</returns>
        public double Sample()
        {
            return DoSample(RandomNumberGenerator, mMean, mStdDev);
        }

        /// <summary>
        /// Samples an array of normal distributed random variables.
        /// </summary>
        /// <param name="size">The number of variables needed.</param>
        /// <returns>An array of random numbers from this distribution.</returns>
        public double[] Sample(int size)
        {
            return DoSample(RandomNumberGenerator, size, mMean, mStdDev);
        }

        #endregion

        /// <summary>
        /// Checks the parameters of a normal distribution.
        /// </summary>
        /// <param name="stddev">The standard deviation of the normal distribution.</param>
        /// <exception cref="ArgumentOutOfRangeException">If the standard deviation is negative.</exception>
        private static void CheckParameters(double stddev)
        {
            if (stddev < 0.0)
            {
                throw new ArgumentOutOfRangeException("stddev", Resources.StandardDeviationCannotBeNegative);
            }
        }

        /// <summary>
        /// Evaluates the cumulative distribution function for the normal distribution.
        /// </summary>
        public double CumulativeDistribution(double x)
        {
            return 0.5*(1.0 + SpecialFunctions.Erf((x - mMean)/(mStdDev*(System.Math.Sqrt(2.0)))));
        }

        /// <summary>
        /// Evaluates the inverse cumulative distribution function for the normal distribution.
        /// </summary>
        public double InverseCumulativeDistribution(double p)
        {
            return mMean - mStdDev*System.Math.Sqrt(2.0)*SpecialFunctions.ErfcInv(2.0*p);
        }

        /// <summary>
        /// Samples a pair of standard normal distributed random variables using the Box Muller transform.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="r2">The second random number.</param>
        internal static double SampleBoxMuller(System.Random rnd, out double r2)
        {
            double v1 = 2.0*rnd.NextDouble() - 1.0;
            double v2 = 2.0*rnd.NextDouble() - 1.0;
            double r = v1*v1 + v2*v2;
            while (r >= 1.0 || r == 0.0)
            {
                v1 = 2.0*rnd.NextDouble() - 1.0;
                v2 = 2.0*rnd.NextDouble() - 1.0;
                r = v1*v1 + v2*v2;
            }
            double fac = System.Math.Sqrt(-2.0*System.Math.Log(r)/r);
            r2 = v2*fac;
            return v1*fac;
        }

        /// <summary>
        /// Samples a normal distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="mu">The mean of the normal distribution.</param>
        /// <param name="sdev">The standard deviation of the normal distribution.</param>
        public static double Sample(System.Random rnd, double mu, double sdev)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(sdev);
            }
            
            return DoSample(rnd, mu, sdev);
        }

        /// <summary>
        /// Samples an array of normal distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of variables needed.</param>
        /// <param name="mu">The mean of the normal distribution.</param>
        /// <param name="sdev">The standard deviation of the normal distribution.</param>
        public static double[] Sample(System.Random rnd, int n, double mu, double sdev)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(sdev);
            }

            return DoSample(rnd, n, mu, sdev);
        }

        private static double DoSample(System.Random rnd, double mu, double sdev)
        {
            double v2;
            return mu + sdev*SampleBoxMuller(rnd, out v2);
        }

        private static double[] DoSample(System.Random rnd, int n, double mu, double sdev)
        {
            double[] arr = new double[n];
            for (int i = 0; i < n; i++)
            {
                arr[i] = DoSample(rnd, mu, sdev);
            }
            return arr;
        }
    }
}